setwd('/Users/emmatimminsschiffman/Documents/Dissertation/secondary stress/glycogen')
glyc<-read.csv('glycogen.csv', header=T)

summarySE <- function(data=NULL, measurevar, groupvars=NULL, na.rm=FALSE,
                      conf.interval=.95, .drop=TRUE) {
    require(plyr)

    # New version of length which can handle NA's: if na.rm==T, don't count them
    length2 <- function (x, na.rm=FALSE) {
        if (na.rm) sum(!is.na(x))
        else       length(x)
    }

    # This does the summary. For each group's data frame, return a vector with
    # N, mean, and sd
    datac <- ddply(data, groupvars, .drop=.drop,
      .fun = function(xx, col) {
        c(N    = length2(xx[[col]], na.rm=na.rm),
          mean = mean   (xx[[col]], na.rm=na.rm),
          sd   = sd     (xx[[col]], na.rm=na.rm)
        )
      },
      measurevar
    )
# Rename the "mean" column    
    datac <- rename(datac, c("mean" = measurevar))

    datac$se <- datac$sd / sqrt(datac$N)  # Calculate standard error of the mean

    # Confidence interval multiplier for standard error
    # Calculate t-statistic for confidence interval: 
    # e.g., if conf.interval is .95, use .975 (above/below), and use df=N-1
    ciMult <- qt(conf.interval/2 + .5, datac$N-1)
    datac$ci <- datac$se * ciMult

    return(datac)
}

plot.glyc<-summarySE(glyc, measurevar='Glycogen', groupvars='Treatment')
#    pCO2 N Glycogen        sd        se        ci
#1  400 8 1.444330 0.8772902 0.3101689 0.7334329
#2  800 8 1.122244 0.5144093 0.1818712 0.4300570
#3 2800 8 2.139174 1.2277117 0.4340616 1.0263927

pd<-position_dodge(0.1)
ggplot(plot.glyc, aes(x=Treatment, y=Glycogen))+
geom_errorbar(aes(ymin=Glycogen-ci, ymax=Glycogen+ci), width=0.1, position=pd)+
geom_line(position=pd)+
geom_point(position=pd, size=3)+
xlab('pCO2 (µatm)')+
ylab('Concentration Glycogen µg/µl')

glyc.lm<-lm(glyc$Glycogen~glyc$Treatment)
anova(glyc.lm)
#F=2.5507, p=0.1019

#corrected for original tissue mass used for reaction
#each glycogen pellet was reconstituted in 200 µl of water
Glycogen.per.mg<-glyc$Tissue.mass*(200/glyc$Glycogen)
#units is µg glycogen per mg tissue
glyc2<-cbind(glyc, Glycogen.per.mg)
glyc2.lm<-lm(glyc2$Glycogen.per.mg~glyc2$Treatment)
anova(glyc2.lm)
#F=0.9446, p=0.4047

plot.glyc2<-summarySE(glyc2, measurevar='Glycogen.per.mg', groupvars='Treatment')

pd<-position_dodge(0.1)
ggplot(plot.glyc2, aes(x=Treatment, y=Glycogen.per.mg))+
geom_errorbar(aes(ymin=Glycogen.per.mg-ci, ymax=Glycogen.per.mg+ci), width=0.1, position=pd)+
geom_line(position=pd)+
geom_point(position=pd, size=3)+
xlab('pCO2 (µatm)')+
ylab('Concentration Glycogen µg/mg')

Treatment2<-as.numeric(c(rep('400', 8), rep('800', 8), rep('2800', 8)))
glyc2<-cbind(glyc2, Treatment2)
plot(y=glyc2$Glycogen.per.mg, x=glyc2$Treatment2, type='p', xlab='pCO2 (µatm)', ylab='Glycogen µg/mg Tissue', xaxp=c(400, 2800, 6), col=c(rep('blue', 8), rep('forestgreen', 8), rep('orangered', 8)), cex=1.25)

#remove outlier from 2800 
glyc.rm<-(glyc2$Glycogen.per.mg[-20])
treatment.rm<-(glyc2$Treatment2[-20])
lm.rm<-lm(glyc.rm~treatment.rm)
anova(lm.rm)
#p=0.8369
plot(y=glyc.rm, x=treatment.rm, xaxp=c(400, 2800, 6), col=c(rep('blue', 8), rep('forestgreen', 8), rep('orangered', 7)), cex=1.25, ylab='Glycogen content (µg per mg tissue)', xlab='pCO2 (µatm)', xaxp=c(400, 2800, 6))